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Abstract 

This paper deals with the mathematical modelling of large strain 
electro-viscoelastic deformations in electro-active polymers. Energy 
dissipation is assumed to occur due to mechanical viscoelasticity of 
the polymer as well as due to time-dependent effective polarisation of 
the material. Additive decomposition of the electric field E = Eg -|- 
and multiplicative decomposition of the deformation gradient F = 

FgF^ are proposed to model the internal dissipation mechanisms. The 
theory is illustrated with some numerical examples in the end. 
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1 Introduction 

Over the past few decades, the theory and nnmerical simnlation of the con- 
pled electro-mechanical problem have been interesting snbjects of research, 
cf. Pao [1] and Eringen and Mangin |2]. However, with the invention of the 
so-called electro-active polymers (EAPs) capable of exhibiting large deforma¬ 
tions in response to the application of electric helds, several new challenges 
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appear and need to be addressed. Open problems remain both in the nnder- 
standing of electro-mechanical conpling in soft matter and in simnlating the 
behavior of electro-sensitive bodies nnder the inflnence of an electric held. 

EAPs can be nsed as alternatives to materials traditionally nsed to de¬ 
velop actnators like piezoelectric ceramics, shape memory metals and electro- 
rheological hnids, cf. O’Halloran et al. |3]. Potential applications of EAPs 
in developing artihcial mnscles and robotic systems inclnde robot manipn- 
lators |1], soft pnmps j5], lond-speakers [6], portable force feed-back devices 
[7], haptic interfaces [8], electric generators for energy harvesting [^-|10]. 
transport vehicles mi. mi. and sensing eqnipment [I3]^[I1], among others. 

Efforts were made in the past to model and simnlate the behavionr of 
EAPs nsing the theory of nonlinear elasticity and nonlinear visco-elasticity, 
for example, by Kofod [15], Sommer-Larsen et al. [T6|, Gonlbonrne et al. mi. 
Yang et al. [HI |T9|, and Rosset et al. [20]. However, the papers mentioned 
above assnme that the material electric properties are independent of defor¬ 
mation. Note that becanse large strain occnrs dnring the deformation pro¬ 
cess, the nonlinearity of the material electric properties must be accounted 
for. In order to overcome this shortcoming, some related boundary-value 
problems involving hnite deformation were analyzed by taking into account 
the nonlinearity of the electric polarisation, for example, in the works of 
Voltairas et al. [21], Dorfmann and Ogden [22], Muller et al. [22], Zwecker 
et al. [21], and Vertechy et al. [22]. The effect of viscosity in the modeling of 
EAPs was examined recently by Ask et al. [26] and Biischel et al. mi. 

A basic assumption in the modelling of EAPs in the papers mentioned 
above has been an instantaneous or ‘elastic’ response of the material to an 
applied electric held. This, however, may not be the case in all the electroac¬ 
tive polymers and we aim at modelling this phenomena in this research. We 
work under a more general case where it is assumed that on the application 
of an electric held, the overall macroscopic polarisation of the material is 
time-dependent. Thus, in addition to the mechanical viscoelasticity of the 
polymeric matrix, an additional energy-dissipating mechanism is considered 
on account of the evolution of the electric polarisation with time. 

Among the several approaches towards a phenomenological theory of me¬ 
chanical viscoelasticity, the literature is usually divided on account of the 
nature of internal variable used to quantify dissipation. The internal vari¬ 
able can be assumed to be of stress-type, as proposed by Simo [21] and Lion 
[22], or it can be strain-type, as used by Lubliner [SU], Reese and Govind- 
jee EH and Huber and Tsakmakis [32]. In the latter approach, which has 
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been also followed in this paper, the deformation gradient is decomposed 
into elastic and inelastic parts (F = FeF^) where the inelastic part is de¬ 
termined from a differential type flow rule. In addition to the mechanical 
viscoelastic dissipation, we also model electric dissipation by considering a 
similar decomposition of our independent variable (electric field in this case) 
into ‘elastic’ and ‘viscous’ parts as E = Eg -|- E^. This follows a similar ap¬ 
proach by Saxena et ah [33] for the magnetic counterpart of this problem. 
The energy and momentum balance laws of electroelasticity are derived from 
the fundamental equations of electrostatics following the work of McMeeking 
and Landis [M] . 

This paper is organised as follows. The theory of rate-dependent elec¬ 
troelastic deformations is presented in Section 2. Starting with the basic 
principles of electrostatics and continuum mechanics (as detailed in the Ap¬ 
pendix), we obtain the energy and momentum balance laws in the case of 
electroelasticity. The deformation gradient and the electric field are decom¬ 
posed into equilibrium and non-equilibrium parts (F = FeF.i;,E = Eg -|- E^). 
Using the laws of thermodynamics and a form of the free energy density 
function, constitutive equations are derived along with the conditions to be 
satisfied by the evolution equations of the non-equilibrium quantities. 

For the purpose of obtaining numerical solutions later, the energy density 
function and the evolution equations for the non-equilbrium quantities are 
specialised to specific forms. Several electro-visco-elastic coupling parameters 
are introduced in this step and we define thermodynamically consistent and 
physically reasonable evolution laws for the internal variables. In Section 
3, numerical solutions are obtained corresponding to five different types of 
(mechanical and electric) loading conditions. The effects of the underlying 
deformation, strain rate, electric field, and electric field rate are studied on 
the evolution of the resulting stress and the dielectric displacement. The 
results, presented graphically, show a strong coupling between strain and 
electric field, as well as the strong dependence of the response on electro¬ 
viscoelastic coupling parameters thus making the model amenable to fitting 
with experimental data, as an when it becomes available in future. 


2 Theory 

We consider an electroelastic material that, when undeformed, unstressed 
and in the absence of electric fields, occupies the material configuration Bq 
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with boundary BEq. It is then subjected to a static deformation due to the 
combined action of an electric held, mechanical surface tractions and body 
forces. The spatial conhguration at time t is denoted by Et with a boundary 
dEt- The two conhgurations are related by a deformation function x which 
maps every point X G to a point x = x(X, t) G Ef. The deformation 
gradient is dehned as F = Gradx, where Grad is the gradient operator with 
respect to X. Its determinant is given by J = det F. 

2.1 Balance laws and boundary conditions 

2.1.1 Equations of electrostatics 

Let q be the electric charge density per unit volume in Bt, e be the spatial 
electric held vector, d be the spatial electric displacement vector, and p 
be the spatial polarisation vector. The balance equations for the electric 
quantities are given by a simplihed form of the two Maxwell’s equations as 

curie = 0, divd = g, (1) 

where the electric vectors are related by the constitutive law 

d = £o® + P, (2) 

and curl and div denote the corresponding diherentiation with respect to the 
position vectors x in the spatial conhguration Et. 

We note that the above equations can also be written in the material 
conhguration Bq by employing the following transformations 

E = F'e, D = JF^M, P = JF'^p, (3) 

thus giving 

GurlE = 0, DivD = Jg, D = + P, (4) 

such that Gurl and Div denote the corresponding diherentiation operators 
with respect to the position vectors X in Bq and C = F^F is the right 
Gauchy-Green deformation tensor. 

Since curl of a gradient vanishes, the electric held vector can be written 
as the gradient of a scalar potential from equation (j^i as 

(E = —grad (j). (5) 
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At an interface or a boundary, the electric vectors must satisfy the con¬ 
ditions 

n X |e] = 0, n ■ |(d] = g, (6) 

where q is the surface charge density, n is the unit outward normal to the 
surface and !•] represents the difference . 


2.1.2 Linear and angnlar momentum balance 

The balance of linear momentum in the conhguration Bt is given in terms of 
the total Cauchy stress tensor as 

div cr*°*fm = pa. (7) 

Here is the total Cauchy stress tensor that takes both mechanical and 
electric effects into account, fm is the purely mechanical body force, p is the 
mass density, a is the acceleration, and the divergence operator is taken to 
operate on the hrst index of a second order tensor. We refer to Appendix A 
for a detailed derivation of the balance equations in the context of electroe¬ 
lasticity. 

The above equation can be written in referential form using the total 
Piola-Kirchhoff stress as 


Div (S*°*F*) + fM = p.a. 


( 8 ) 


with pr = Jp being the referential mass density and Im = Jim being the 
referential body force. Note that the tensor S = JF^^crF'* is sometimes 
also referred to as the ‘second’ Piola-Kirchhoff stress. 

The principle of balance of angular momentum renders the Cauchy and 
the Piola-Kirchhoff stress tensors symmetric 


CT“‘)‘ = cr“, (S“)‘ = S“*. 


(9) 


The corresponding boundary conditions are given by the equations (53) 


and (54). 


2.1.3 Internal dissipation 

Very often, the EAPs are synthesized from a rubber like polymer. The poly¬ 
meric rubber matrix is viscoelastic in nature which leads to energy dissipation 
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on a mechanical deformation. In addition to this, energy dissipation can also 
occur due to a time-dependent polarisation of the material on application of 
an electric held. We consider the possibility that on a sudden application of 
an electric held (or a potential diherence), the electric displacement D, the 
polarisation P, and the resulting electric contribution to stress generated in 
the material evolve with time to reach an equilibrium value. Thus, the two 
ehects need to be modelled appropriately. 

To take into account mechanical viscous ehects, we assume the existence 
of an intermediate conhguration Bi that is, in general, incompatible. The 
tangent spaces of Bq and Bi are related by a second order tensor that 
quantihes viscous motion while the tangent spaces of Bi and Bt are related by 
a second order tensor Fe that quantihes elastic distortion. The conhguration 
Bi is in parallel to the energy-conserving electroelastic deformation from Bq 
to Bt- This motivates the decomposition of the deformation gradient into an 
elastic and a viscous part (cf. Lubliner [30] and Reese and Govindjee [3T] ) 
as 

F = FeF,. (10) 

For future use we dehne the right and the left Cauchy-Green strain tensors 
as C = F*F and b = FF*, respectively. Corresponding quantities are dehned 
for the viscous parts as Cy = F*F^ and = F^F*. 

In order to model the electric dissipation ehects, we consider an additive 
decomposition of the electric held vector into an ‘elastic’ and a ‘viscous’ part 

<E = (Be -|- E = Eg “1“ E^, (H) 

which dehnes a dissipation mechanism in parallel to an energy-conserving 
mechanism that gives the instantaneous response. The above additive de¬ 
composition of the electric held is motivated by a similar decoupling of the 
deformation gradient into elastic and viscous parts in the viscoelasticity the¬ 
ory. A similar approach has been taken by Saxena et ah [33] to model 
magneto-viscoelastic ehects wherein an additive decomposition of the mag¬ 
netic induction as B = Bg -|- has been performed. The behaviour of the 
internal variables dehned above is assumed such that if a constant electric 
held is applied at time t = 0, then at that instant E^ = 0 and Eg = E. As 
time progresses, E„ —)■ E and Eg —)■ 0. 

The material response is assumed to be given by an energy density func¬ 
tion G that depends on the parameters C, C„,E,E^ and the temperature 
(see Appendix A for details). Taking partial derivatives of G with respect to 
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its arguments and substituting in the thermodynamical inequality (79), we 
get 




'an 


-h D 

[dE \ 
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m + 
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on 


^Q-Grad79>0, (12) 


where p is a Lagrange multiplier associated with the constraint of incompress¬ 
ibility and for the sake of brevity we now refer to the total Piola-Kirchhoff 
stress as simply S. A superposed dot represents a material time deriva¬ 
tive. 

On an application of the procedure due to Coleman and Noll |35], we 
arrive at the following constitutive relations for the total Piola-Kirchhoff 
stress, the electric displacement and the specihc entropy 


S 




and the dissipation conditions 


dQ 1 dQ 

5E ’ ^ pr d'd' 


dVL • dVL • 
dc,'' " ^ aE, ■ 


^Q-Grad^9<0. 

V 


(13) 


(14) 


It is noted here that if the incompressibility constraint (J = 1) is not 
imposed, then the constitutive equation for stress is simply 


S 



(15) 


It is further assumed that the ‘viscous’ electric held, viscous strain and 
the temperature are independent of each other, thereby reducing the above 
inequality to the following separate conditions 


dQ 


: C, < 0, 


dn 

cIE,, 


■ E„ < 0, 


Q ■ Grad'd < 0. 


(16) 


These conditions must be satished by all processes to be thermodynamically 
admissible. Isothermal conditions are assumed henceforth in this paper and 
we remove the dependence on the temperature 

A common technique useful for performing numerical computation ( j3lj. 
ESI) is to decompose the total energy into an equilibrium part associated 
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with the direct deformation from Bo to Bt, and a non-equilibrium part due 
to the internal variable and the elastic deformation from Bi to Bt- 


fl (C, C^, E, E^,) = fie (C, E) -|- (C, C^, E, E^,). 


(17) 


The above simplihcation is used to write the dissipation inequalities (16) 1,2 


as 


dfly 


: C, < 0, 


dfl^ 

dE„ 


■ E„ < 0. 


(18) 


2.2 Specialised constitutive laws 

To obtain numerical results in order to understand the physical implications 
of the several aspects of the developed theory, we use a prototype energy 
function that is a generalisation of the isotropic neo-Hookean function to 
electroelasticity and is given by 

fie = ^ [Jl — 3] -|- me /4 -|- Ug/s, (19) 

where the scalar invariants Ji, I 4 and J 5 are dehned as 

Jl = C : I, /4 = [E 0 E] : I, Jg = [E O [C’^E] ] : I. (20) 


The parameter pe is the shear modulus in the absence of electric helds, rUe 
and He are electroelastic coupling parameters with rrie/eo and Ue/so being 
dimensionless, and I is the identity tensor in Bo- 

For the non-equilibrium part of the stored energy, a similar (but with 


one different electro-mechanical coupling term) form as that of fie in (19) is 
chosen in a way that the energy depends only on the non-equilibrium ‘elastic’ 
variables Ce = and Ee = E — E^. The reason and the procedure of 

such a simplihcation is discussed in 


fl, = ^[C-. - 3] + m4 [E - E„] 0 [E - E,] ] : I 


+ 


[C [E - E„] ] 0 [C [E - EJ ] 


: I. 


( 21 ) 


Here /r^, m^ and riy are the non-equilibrium parameters with dimensions sim¬ 
ilar to their electroelastic counterparts dehned in (19). 









The evolution equation for the ‘viscous’ electric held vector is taken to be 


=-TT^, 

= [niyl + [E - E^,], 

4-^e 

where we have taken Tg as a time constant to measure evolution of ‘viscous’ 
electric held and is a scaling parameter. For the viscous deformation, we 
use the evolution equation given by Koprowski-Theiss et ah [SZj that is a 
simplihcation of a general form used by Lion |29] 


( 22 ) 

(23) 


c„ = L 

T, 


C - - [C : Cy] C. 


(24) 


Both the evolution equations (23) and (24) satisfy the thermodynamic con¬ 
straints (18 )i^ 2- They are also dehned such that the viscous quantities E^ 
and stop evolving as the system reaches equilibrium. 

For the above dehned energy density functions, the total Cauchy stress 
<T = J-iFSF* is given as 


cr = + cr'’+ pi, (25) 

where i is the identity tensor in Bt and 

cr® = /ieb — 2ne<B ® e, (26) 

cr" = + 2ny [[b ©e] ® [b^eg] + [b^Ce] ® [b ©e]] • (27) 

The total electric displacement d is given as 

d = d" + d^ (28) 

d® = —2meb© — 2ne©, (29) 

d'^’ = —2m^b©e — 2nyh^ee. (30) 


We now use the above dehned energies and specialised constitutive laws 
to perform numerial calculations in the following sections. 
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3 Numerical examples 

In order to understand the physical behaviour predicted by the model de¬ 
veloped in the sections above, we consider some numerical examples. By 
prescribing the deformation and the electric held at a point, we calculate 
the evolution of stress and electric displacement in the material. The depen¬ 
dence of the total Cauchy stress and the electric displacement on modelling 
parameters, applied stretch, and electric potential difference or electric held 
is demonstrated graphically. 

3.1 No deformation 

In the hrst case, in order to highlight the ehects of non-equilibrium electric 
held, we consider no deformation but existence of an electric held. Let the 
test specimen be held hxed (C = I) using appropriate boundary tractions 
and a sudden but constant electric held is applied at time t = 0. This leads 
to the generation of a viscous overstress and an electric displacement, both 
of which settle down to equilibrium values with time. The following values 
of the material parameters are used for calculations in this case 


/ie = 5 X 10® MPa, = 2 X 10® MPa, rrie = —10 N/V^, ^ = 1 

Tie = —6 N/V^, Te = 1 S, Ty = 50 S, = — rUe, Uy = —Ue- (31) 


The electric held is given by 


IE2 — IE3 


0 , for f < 0, 

3 X 10^ V/m, for t>0. 


(32) 


Time integrations of the evolution equations are performed using a modi- 
hed Runge-Kutta scheme implemented in the ode45 solver of MATLAB. It is 
observed that an ‘electro-viscous’ stress is developed in the material at time 
t = 0 when the electric held is switched on. The stress decays with time 
due to evolution of the ‘viscous’ electric held and reaches an equilibrium 
value. The rate of decay and the initial non-equilibrium value depend on 
the material parameters rriy and Uy. As observed from Fig. [^c), changing 
the value of ruy does not change the initial value of stress an but a higher 
TUy causes a faster decay of stress to its equilibrium value. It’s seen from 
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Fig. 1. Evolution of (a,b) electric displacement di (Nm“^V“^) and (c,d) the 
total Cauchy stress an (N/m^) with time for two different values of and 
n„: (a,c) (i) rriv = 3 N/V^, (ii) my = 10 N/V^; (b,d) (i) Uy = 3 N/V^, 
(ii) Uy = 10 N/V^. 
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[^d) that a higher value of results in a higher value of initial stress while 
slowing the decay process. 

Similar to the non-equilibrium stress, a ‘viscous’ electric displacement is 
developed in the material at f = 0 which evolves with time to achieve an 
equilibrium value. As is seen from Figs. Sa) and [^b), a higher value of 
either or leads to a smaller initial value of electric displacement. A 
higher rriv leads to a faster evolution while the reverse happens for a higher 

Uy. 

3.2 Uniaxial deformation 

In this case, we apply both the uniaxial tension and the electric held at f = 0 
as a step function. The stretch Ai and the electric held Ei are both applied 
in the same direction. 

The viscous overstress decays at two time scales - the mechanical part 
corresponding to evolution of C„ and the electric part corresponding to evo¬ 
lution of as shown in Fig. |^a) and[^b) where the curves are plotted for 
two diherent values of the initial deformation Ai = 1.5 and Ai = 1.6, respec¬ 
tively. As expected, a higher Ai results in a higher value of stress but it also 
results in a faster evolution of an. As shown, the end point of a curve in 
Fig-i a) has the same value as the starting point of the corresponding curve 
in Fig. [^b) since they both correspond to the same loading conditions. 

3.3 Time dependent electric field 

In this section, we study the ehect of a time-varying electric held on the 
induced total stress and the overall dielectric displacement in an undeformed 
material. The specimen is taken to be hxed at zero deformation and an 
electric held is applied at time t = 0 in the xi direction with a constant 
rate until a value ei = 300 V/m is obtained. The electric held is then 
reduced with the same rate until it reduces back to zero. Numerical results 
for this case are shown in Figs. and and are obtained for the values 
<Bi = {±200, ±400} V/m-s in Fig. |^and <Bi = ±300 V/m-s in Fig. 

It is observed from Fig. [^a) that the electric displacement increases with 
an increasing electric held in a nonlinear fashion and returns back through 
a diherent path, eventually reaching a non-zero value of d for zero e. If the 
electric held is held at zero after this time, the electric displacement would 
gradually relax to zero showing a behaviour similar to that in Fig. ia). A 
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Fig. 2. Uniaxial loading: Evolution of the total Cauchy stress an (N/m^) 
with time t (s) for (a) short time-scale and (b) large time-scale, (i) Ai = 1.5, 
(ii) Ai = 1.6. 


higher rate (<Bi) of the electric held leads to lower extreme values of d but a 
larger area enclosed within the curve leading to a higher energy loss per cycle. 
Another key observation is that the maximum value of dielectric displacement 
d is obtained on the unloading curve just at the beginning of the unloading 
part of the cycle. This is because the evolution of slightly lags behind 
the instantaneous changes in e. The total Cauchy stress in Fig. |^b) varies 
almost quadratically with e, slight hysteresis being observed only at higher 
rates of electric held. It is also observed that the value of an on the loading 
curve is usually higher than that on the unloading curve for larger values of 
© 1 . 

Similar curves for diherent values of the material parameters riy and 
are plotted in Figs. |^a,b) and|^c,d), respectively. 

A higher Uy increases the maximum value obtained by d and results 
in a larger area under the curve thus increasing the energy dissipated per 
cycle. The Cauchy stress an remains largely unchanged with the value of Uy. 
Increasing the value of lowers the peak value of the dielectric displacement 
but increases that of the total stress. A higher increases the amount of 
hysteresis in each case. 
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(b) 


Fig. 3. Variation of (a) the electric displacement di and (b) total Cauchy 
stress (Til (N/m^) with electric held ©i: (i) <Bi = ±200 V/m-s, (ii) <Bi = 
±400 V/m-s. 


3.4 Time dependent shear 

In this case, we consider the deformation to be given by 

r = R, 9 = Q + aZ, z = Z, (33) 


This is an effective deformation obtained by torsion of a cylinder by constrain¬ 
ing the plane cross-sections to remain plane. It is typically obtained exper¬ 
imentally in a disc deformed using a rheometer. Locally, the deformation is 
a plane strain where a depends on time and is given by o; = aosincat. The 
applied Lagrangian electric held is in the z direction given as E = {0, 0, E^y. 
Let (1,2,3) correspond to {R,Q,Z) in the cylindrical coordinate system, 
then the deformation gradient F and the left Cauchy-green strain tensor b 
are given as 


■ 1 

0 

0 ■ 


■ 1 

0 

0 ■ 

0 

1 

7 

, [b] = 

0 

1±7" 

7 

0 

0 

1 


0 

7 

1 


(34) 


Here 7 = aR is the amount of shear at the radius R. 


Substituting these values in the expressions (25)-(27) for the total Cauchy 
stress, we obtain 


cr23 = l^el + - Iv] - ‘irie'jE^ ± 2n^7^[-3 ± 37 ^ ± 27 '^]F;/ 3 , (35) 
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(c) (d) 

Fig. 4. Variation of (a,c) the electric displacement di (Nm”^V“^) and (b,d) 
the total Cauchy stress an (N/m^) with electric held: (a.b) (i) = 1 N/V^, 

(ii) Uv = 20 N/V^; (c,d) (i) Tg = 5 s, (ii) Tg = 10 s. 
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(a) (b) 


Fig. 5. Shear Cauchy stress CT 23 (N/m^) vs the shear strain 7 . (a) Switching 
on electric held during a dynamic shear experiment. The arrow shows the 
point where electric held is switched on resulting in a jump in (T 23 . (i) E = 
0, (ii) E 3 = 5 X 10^ V/m. (b) Dependence on the electroelastic coupling 
parameter Ue- (i) rie = —1 N/V^, (ii) Ug = —6 N/V^, (iii) Ug = —10 N/V^. 


where 7.0 is the corresponding shear strain for the viscous part of the defor¬ 
mation and = E 3 — Ey^. The evolution equation (24) is reduced to 


d'^y 

dt 


1 

^ V 



[3 + It - 7„]^] T„ 


(36) 


Numerical results are obtained for following values of the material pa¬ 
rameters 


T^ = ls, a; = 2s"\ Tg = 2s, E^ = b x /m. (37) 

In the dynamic shear experiment, when measured over time, stress-strain 
curves form a hysteresis loop (also called the Lissajous plots) with the area 
enclosed by the curve related to the total energy loss per cycle. We show 
the dependence of such a curve on an applied electric held in Fig. |^a). The 
curves in this hgure correspond to switching on of an electric held while the 
dynamic shear test is being performed - curve (i) corresponds to e = 0 while 
curve (ii) corresponds to e = 500 V/m. As observed, the Cauchy stress (< 723 ) 
gets a jump in its value on a sudden application of electric held and then 
evolves due to dissipation to converge to a steady-state region. The ehect of 
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Fig. 6. Dependence of the Lissajous plots on (a) the parameter T^: (i) = 

0.25 s, (ii) Ty = 0.5 s, (hi) = 1 s; (b) the parameter (i) = 2 x 10® 

MPa, (ii) /i„ = 4 X 10® MPa, (hi) /Xy = 6 x 10® MPa and (c) the oscillation 
frequency uj\ (i) a; = 1 s“^, (ii) oj = 2 s“^, (hi) a; = 3 s“^. 

the electroelastic coupling parameter Ue on the dynamic shear experiment is 
demonstrated in Fig. [^b). Increasing the magnitude of Ue causes a stronger 
coupling with the electric held, thus increasing the slope of the curve which 
represents an increase in the effective shear modulus. 

Curves for different values of the parameter are shown in Fig. [^a). 
Increasing the value of results in lowering the slope and hence the effective 
shear modulus. It also results in a larger area enclosed in the curve thus 
increasing the energy loss during each load cycle. Reverse happens in the 
case when the parameter /iy is changed in Fig. |^b). Increasing the value of 
/iy results in a higher slope curve but also more dissipation per load cycle. 
Plot for different oscillation frequencies u in the presence of an electric held 
is shown in Fig. [^c). Increasing the oscillation frequency increases the slope 
of the ellipse thereby increasing the effective shear modulus. 

Dynamic shear test performed in a rheometer is a very widely used ex¬ 
perimental procedure to characterise viscoelastic materials. The curves pre¬ 
sented in Figs. and show that the presented model is able to capture 
serveral aspects of the Lissajous curves and therefore should be able to ht 
the experimental data effectively, as and when it becomes available. 
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3.5 Extension and inflation of a compressible hollow 
cylinder 

Let an infinitely long hollow cylinder of internal radius A and external radius 
B be inflated and stretched in the axial direction such that the new internal 
and external radii are a and b, respectively. Stretch in the axial direction is 
uniform such that a point at an axial position Z is now at the axial coordinate 
z = XzZ in the deformed conhguration. It is assumed that there is no 
variation of deformation with the 0 coordinate and thus the deformation and 
the deformation gradient tensor in the cylindrical polar coordinate system 
{R, 0, Z) are given as 


r{R) = R + u{R), 9 = Q, z = X,Z, 


(38) 


1F| 


g 0 0 

r 

« s « 

0 0 


(39) 


with g = dr/dR. 

Along with the imposed mechanical boundary conditions, a potential dif¬ 
ference A(j) = (j)b — (j)a is applied across the radial direction. This results in 
the generation of an electric held that is coupled with the underlying de¬ 
formation. In order to obtain the exact deformation and the variation of 
the electric held along the radial direction, the following system of equations 
needs to be solved: 


DivD = 0, Div (SF*) =0, E = -Grad0, 
r{A) = a, r{B) = b, 0(A) = 0^, 0(5) = <f>b- 


(40) 

(41) 


We consider a compressible material for this problem and thus the energy 
function (19) is generalised to 

Ge = y [/l - 3] - Pe In J -F y [lu jf + TUeh + Ueh, (42) 


while the energy for the viscous part is still given by (21). 

Due to the geometric symmetry of the cylinder and the existence of po¬ 
tential diherence across the radial direction, the electric held vanishes in the 
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axial and the circumferential directions. For the given energy function, the 
electric displacement in radial direction is given as 


Br = -2 [irieER + Ueg + myEe{R) + Uyg'^EeiR)] , (43) 

while the nominal stress T = SF* in the three principal directions is given 
as 

1 I! (1 ^T) 

Trt = g^eg + [~ge + Ae 111 J] + —- ^ + 2nyE^^^^~^g^, (44) 

2g ^v{RR) g 

Tee = ge-j^ + ^ [—ge + Ae In J] + -, (45) 

H 2r KUy(^ee) 

Tzz = ge^z + -jr [~ge + Ae lu J] + gy— -. (46) 

^ ^v{ZZ) 


The three principal components in the radial, azimuthal and axial directions 
of the tensor are denoted by Cy(RR),Cy(ee) and Cy(zz), respectively, in 
the above equations. 

The governing equations ([40|) are reduced to 


a(i?Dii) 

dR 


djRTRR) 

dR 


T, 


00 , 



(47) 


which need to be solved along with the boundary conditions (41) for every 
time step. The values of Eg and are updated at every time step using 
the solution for r{R),(j){R), and the evolution equations (23) and (24). It is 
noted that none of the governing equations relevant to this problem depend 
on the axial stretch A 2 . 


3.5.1 Numerial calculations 

Numerical calculations are performed for the following values of the material 
parameters and physical quantities 

/le = 5 X 10®N/m^ Ae = 6.6667 X 10® N/m^ 
me = —= —10 N/V^, Ue = —0.06N/V^, ny = 6'N/V‘^. (48) 

Computations are performed using a hnite-difference technique implemented 
in the dsolve solver of Maple for solving BVPs. 
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0 5 10 ^ 15 

(a) (b) 

Fig. 7. (a) Displacement AR (mm) along the radius of the tube at time 
t = {0,5,10} s. (b) Pressure (N/m^) at (i) the internal boundary and (ii) 
the external boundary required to maintain the specihed geometry. 





(a) (b) 

Fig. 8. (a) Radial electric held E/j (V/m) vs radius R (mm), (b) Radial elec¬ 
tric displacement Hr (N m“^ V“^) vs radius R (mm) at times t = (0, 5,10} s. 
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Fig. 9. Displacement AR (mm) and the electric field (V/m) along the 
radius of the tube at t = {0, 5,10} s for = 5 N/V^. 


In order to understand the effects of the electric dissipation process, we 
neglect the mechanical viscosity for the time being. At time t = 0, a potential 
difference of A0 = 10 V is applied across the internal and external surface of 
a tube with the internal and external radii given by A = 10 mm and B = 20 
mm, respectively. The resulting values of displacement AR = r — R as a. 
function of radius R, electric field E, and the electric displacement D are 
plotted with time in Figs. 00 

In the first case, the internal and external boundaries of the tube are 
fixed such that a = A and b = B. The curves in Fig. [^a) correspond to 
the displacement of points at different instants of time. Due to evolution of 
the ‘viscous’ electric field E^ the displacement changes, starting from a small 
initial value it attains a larger equilibrium value with time. In Fig. [^b), we 
plot the evolution with time of the internal and external pressures required 
to maintain the specified boundary conditions. 

For the same problem, variation of electric field along the radius of the 
tube is plotted in Fig. |^a) at different values of time. It is observed that the 
spatial gradient of the electric field decreases with time. There also exists 
a value of radius (say R = Rq) at which the value of electric field remains 
constant and for R > Rq the magnitude of electric field increases while for 
R < Ro the magnitude of electric field decreases with time. The electric 
displacement attains maximum value at the internal boundary and mini- 
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Fig. 10. Displacement AR (mm) vs radius R for at t = {0,5,10} s for 
(a) b = 19.5 mm, (b) & = 19 mm. 


mum value at the external boundary. Evolution of causes D to increase 
uniformly with time. 

The effect of electro-viscoelastic coupling parameter is shown in Fig. 

in which an smaller value of = 5 Pa m^/V^ is taken while keeping other 
parameters and boundary conditions the same. The resulting mechanical 
displacements and electric fields are plotted in Figs, [^a) and [^b) , respec¬ 
tively. It is observed in this case that the peak value of deformation is higher 
than for = 10, and also this value is approached faster in comparison to 
the previous case. Similar observations are made corresponding to the elec¬ 
tric held. The peak values attained are higher in magnitude with a higher 
gradient along the radius. 

As a last illustration, the mechanical boundary conditions for the prob¬ 
lem are changed keeping other conditions constant. The internal wall of the 
cylinder is held at the same position while the external is compressed and 
brought to a known hxed radius. We consider the displacement of the exter¬ 
nal wall by 0.5 mm and 1 mm (corresponding to 6 = 19.5 mm and 6 = 19 
mm, respectively) and plot the displacement along the radius in Figs. 10 (a) 
and mb) respectively. In this case, as the prescribed mechanical deforma¬ 
tion of the outer wall is increased, the overall variation of AR with R tends 
to be linear. Moreover, the evolution of E^ also tends to have a lower effect 
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on changing the values of AR in the case of larger values of prescribed dis¬ 
placements. This happens since in this case mechanical effects dominate the 
electrostatic forces. 

Acknowledgements: This work is supported by an ERC advanced inves¬ 
tigator grant towards the project MOCOPOLY. 


A Fundamental equations 

In the following sections, we derive the governing equations and the con¬ 
stitutive laws of electroelasticity from the basic principles of electrostatics 
following the approach of Mcmeeking and Landis |34j . 


A.l Balance laws 

A. 1.1 Linear and angular momentum balance 

Conservation of mass implies that for any volume V, the following relation 
must hold 

J p dV = 0 p + p div V = 0, (49) 

V' 

where p(x, t) is the mass density of the material and v is the velocity of a 
point X at time t. A superposed dot implies the material time derivative. 
The laws of balance of linear and angular momenta are given as 

J [im + ffe] dY -h y [tm + tg] dS = j pw dV, (50) 

V S V 


j X X [fm -f- ffe] dv + 


X X [tm + te] dS 


y X X p'vdV. 


(51) 


V S V 

Here, and ffe are, respectively, the mechanical and the electrical body 
force per unit mass, and tm and te are, respectively, the mechanical and 
the electrical surface tractions. The electric body force is considered to be a 
result of the electric held acting inside the material and can be derived from 
a Maxwell (or electrical) stress tensor cTg such that 


ffe = div CTe. 


(52) 
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Thus, the electrical traction on the surface S is given by 

= [cTeJn, (53) 

n being the unit outward normal to the surface S. 

For an inhnitesimally small surface element, in order to balance the linear 
momentum the mechanical Cauchy stress cr^ must balance the mechanical 
and the electrical traction, thus giving 

tm + te - (T^^n = 0, ^ =-{o-ln + o-lJn. (54) 


Using the above equations, the balance of linear and angular momenta 
are now given as 


and 


div {(Tm + (Te) + fm = pa. 

(55) 

CTm + CTe = + ( t \. 

(56) 


A. 1.2 External work on the system 

Rate of work done by the agencies external to the system is given by, cf. [M] 


IV = / f,„-vdU + 


vd^T / 0[gdU]+ / (j) qdS 


(57) 


Bt 


dBt 


Bt 


dBt 


Here the first two terms correspond to the work done by the mechanical 
body force and the surface tractions, respectively. The third term is the 
statement of the rate of electric work - electric potential 0 multiplied by the 
rate of increment of charge qdV. This is the work required in bringing an 
inhnitesimal charge qdV from inhnity to the point where the potential is 
0. The fourth term is the rate of work corresponding to the surface charge 
density. 

It can be shown that 


and 


[qdV] = 


dq 

— + V ■ grad q + q div v 
at 


dU, 


(58) 


[|cl ■ n] dS*] = d + d div v — [grad v] d 


n 


dR. 


(59) 
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In the second equation, use has been made of the Nanson’s formula connect¬ 
ing material and spatial surface elements as ndS* = JF“*N dSr- 


Substituting the above equations in (57) and using the boundary condi¬ 


tion (54), we get 


W= [pa - div (cr^cTe)] ■ V dl/ - / [|cr^-h cr*]n] ■ v dS* 


Bt 


+ 


Bt 


dBt 

dq 

-)- V ■ grad q + q div v 
ot 


dl/ 


+ 


d -I- d [div v] — [grad v] d ■ n dS*. (60) 


dBt 


The two surface integrals can be converted to volume integrals using the 
divergence theorem to give 


W = 


Bt 


[pa - div (cr^ <7^)] ■ v -)- 


dq 

-|- V ■ grad q + q div v 
ot 


div(0 d-|-d [div v] — [grad v] d ) -|-div ([cr^-|-erg] v) did (61) 


On substituting the relations (02, 0, and ( [4^ in the above equation, 
and a subsequent rearrangement of terms, we obtain 


W = 


pa ■ V + [(Tm -|- cTe — d ® e] : grad v -|- e ■ d -|- [div v] e ■ d 


du. 


Bt 


( 62 ) 

We now refer to the expression [cr^ + erg] as the total stress tensor which 


is a symmetric tensor (from equation (56)) including both the mechanical and 
the electric effects. Thus a total symmetric Piola-Kirchhoff stress can now 
be dehned as [cr^ -|- cTg] F~* using which the above statement 

for work rate is written in the referential form as 


W = 


PrV-a + 


tot 


C + E- 


dl4. 


(63) 
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A.2 Thermodynamics and constitutive relations 

By the first law of thermodynamics, balance of energy is stated as a relation 
between work done on the system, internal energy of the system and the heat 
transferred to the system as 

il = W + Q, (64) 


where U represents the energy stored in the system and Q is the rate of heat 
transfer to the system. Specihcally, they can be written as 


d 

dt 


d 

dt 


Bt 


/ T”® ■ 

.7 d n 

•edu+— / -pv-vdu, 
dt J 2 

(65) 

Bt 

Bt 


pr du — 

j q ■ n ds. 

(66) 


Bt 


dBt 


where u is the internal energy per unit mass and r is the rate of heat transfer 
per unit mass. 

In referential form, they give 


if 



PrU + —EoJK ■ [C + —p, 


V ■ V 


dK, 


(67) 


Q= [Prr - DivQ] dK-, 


( 68 ) 


Bo 


where Q = JF-iq is the heat flux in the reference conhguration. 
Thus the statement of hrst law (64) is given in referential form as 


PrU - : C - E ■ P + ^E ■ [eo^C'^E] - ^E ■ L^C'^E 


—prT + Div Q = 0. 


(69) 


We now dehne a total energy function similar to the one given by Dorf- 
mann and Ogden [38] to take into account the energy of electric held 


■0 = PrU — -EoJK ■ [C ^E] . 


(70) 
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On substituting ^|J into the equation above, we obtain 

i) - ; C - E ■ P - + DivQ + [C'^E] ■ E = 0 . (71) 

On using a Legendre’s transformation to change the dependent variable 
from P to E 

^ - E ■ P, (72) 

we eventually obtain the hrst law in the desired form 

(^_ Igtot ; c + D-E-p^r + DivQ = 0. (73) 

Let s be the entropy density per unit mass, then the second law of ther¬ 
modynamics is given as 

(74) 

Bt Bt dBt 

which when written in terms of referential quantities, gives 

j^jprsdV> j ^dK - / (75) 

Br Br dBr 


On using the divergence theorem on the last term, we write the above 
inequality in local form to obtain 


PrS > ^ - ^DivQ -h ■ Grad?9. 
tr V 


(76) 


Substituting the statement of hrst law (73) into the inequality above gives 
1 . 


p + ; C - D ■ E - pr^s - -Q ■ Grad^9 > 0. 

2 V 


(77) 


A hnal Legendre’s transformation is now performed to change the depen¬ 
dency from s to "d 

O = 99 — pr'&s. (78) 


Substitution of the above dehned fl into inequality (77) gives 

— Q -\—: C — D ■ E — prsd —-Q ■ Grad'd > 0. 
2 'd 


(79) 
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We now note that for an incompressible material, the constraint J = 1 
gives 

72 = C-i : C = 0. (80) 

Thus a scalar multiple of this zero term can be added to the above inequality 
without changing its meaning. 
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